%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ESTIMATE THE EFFECT OF TRANSPORTATION COSTS ON DEFORESTATION IN THE AMAZON

% "DEFORESTATION IN THE AMAZON:
% A UNIFIED FRAMEWORK FOR ESTIMATION AND POLICY ANALYSIS"

% by Eduardo Souza-Rodrigues

% This version: November 2018

% OBS: The data set has to be already loaded before running this program.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% THE MODEL:

% ln(Y/(1-Y)) = TC'alpha(U) + X'beta(U) + beta_cons(U)

% where:
% Y     = share of deforested area
% TC    = transportation cost to the nearest port
% X     = covariates
% U     = municipality-level unobservable 
% D     = [D1 D2], straight-line distances to: 
%         (i) nearest port and (ii) nearest capital,
%         D is the IV for TC
% alpha = parameter of interest
% beta  = finite dimensional parameter (covariates)

% Moment Restriction:
% E[ 1{ln(Y/(1-Y)) - TC'alpha(u) - X'beta(u) - beta_cons(u) < 0} - u | D, X] = 0

% Estimation methods: Chernozhukov and Hansen's IVQR

% Unit of observation: municipalitites in the Brazilian Amazon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% (A) PREPARE THE DATA SET

% Name of the variables
Code        = data(:,1);      % Municipality Identifier (IBGE's code)
Y           = data(:,2);      % Proportion of Deforestation
Qa1         = data(:,3);      % Productivity Index - Only Crops
Qa2         = data(:,4);      % Productivity Index - Crops and Pasture
Qrice       = data(:,5);      % Yields Rice
Qman        = data(:,6);      % Yields Manioc (Cassava)
Qcorn       = data(:,7);      % Yields Corn
Qsoy        = data(:,8);      % Yields Soy
Qbean       = data(:,9);      % Yields Beans
TC          = data(:,10);     % Transportation cost to the port
Alt         = data(:,11);     % Altitude
Temp        = data(:,12);     % Temperature
Rain        = data(:,13);     % Rain
Slope       = data(:,14);     % Average Slope
Soil        = data(:,15:18);  % Proportion of soil quality = [prop_potenc2 prop_potenc3 prop_potenc4 prop_potenc5]
Pop         = data(:,19);     % Population
Mining      = data(:,20);     % Presence of Mining
Powerplant  = data(:,21);     % Presence of Powerplant
PPNeighbor  = data(:,22);     % Indicator for Powerplant Neighbor
DistPA      = data(:,23);     % Straight Line Distance to nearest Protected Area
Title       = data(:,24);     % Proportion of the Private Area with Land Title
Fines2005   = data(:,25);     % Fines up to 2005
Fines2003   = data(:,26);     % Fines up to 2003
D           = data(:,27:28);  % Straight Line Distances = [dist to port, dist to state capital]
Ibama       = data(:,29);     % Distance to the Nearest Ibama Agency
CarbonDiff  = data(:,30);     % Difference between Forest Carbon Stock and Deforested Carbon Stock
Area        = data(:,31);     % Area occupied by farms of the size considered
G_IR        = data(:,32);     % Clusters, IBGE Immediate Regions
G50         = data(:,33);     % Clusters, G = 50
G75         = data(:,34);     % Clusters, G = 75
G100        = data(:,35);     % Clusters, G = 100
G125        = data(:,36);     % Clusters, G = 125
G150        = data(:,37);     % Clusters, G = 150
Wpop        = data(:,38);     % Spatially Lagged Population
Wmining     = data(:,39);     % Spatially Lagged Mining
Wpowerplant = data(:,40);     % Spatially Lagged Powerplants
WdistPA     = data(:,41);     % Spatially Lagged Proximity to Protected Areas

% Number of observations
T  = size(Y,1);    

% Truncate Y so that 0 < Y < 1
Yt = Y + (myzero - Y).*(Y<myzero) + (1 - myzero - Y).*(Y>1-myzero);

% Compute log odds ratio
lnY = log(Yt./(1-Yt));

% Spatially Lagged Regressors
WX = [Wpop, Wmining, Wpowerplant, WdistPA];

% Set of Covariates (Model Specification):
if model == 0
    
    % Main Specification
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Title Fines2005 Ibama WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population','Prop Land Title',...
            'Fines','Dist to Ibama','W Population','W Mining','W Powerplant','W Dist to PA','Constant'};
        
elseif model == 1
    
    % No Population
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Title Fines2005 Ibama WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Prop Land Title',...
            'Fines','Dist to Ibama','W Population','W Mining','W Powerplant','W Dist to PA','Constant'};
       
elseif model == 2
    
    % No land title proxy
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Fines2005 Ibama WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population',...
            'Fines','Dist to Ibama','W Population','W Mining','W Powerplant','W Dist to PA','Constant'};
    
elseif model == 3
    
    % No distance to ibama, nor fines
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Title WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population','Prop Land Title',...
            'W Population','W Mining','W Powerplant','W Dist to PA','Constant'};
    
elseif model == 4
    
    % No distance to ibama
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Title Fines2005 WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population','Prop Land Title',...
            'Fines','W Population','W Mining','W Powerplant','W Dist to PA','Constant'};

elseif model == 5
    
    % No Fines
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Title Ibama WX];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population','Prop Land Title',...
            'Dist to Ibama','W Population','W Mining','W Powerplant','W Dist to PA','Constant'};

elseif model == 6
   
    % No Spatially Lagged Regressors
    X     = [Alt Temp Rain Slope Soil DistPA Mining Powerplant PPNeighbor Pop Title Fines2005 Ibama];
    Names = {'TC','Altitude','Temperature','Rain','Slope','Soil 2','Soil 3','Soil 4','Soil 5','Dist to PA',...
            'Mining','Powerplant','Power Plant Neighbor','Population','Prop Land Title',...
            'Fines','Dist to Ibama','Constant'};
    
end

% Dimension of X        
dx = size(X,2);   
   
% All Regressors -- Including TC and the constant term
Xc  = [TC X ones(T,1)];
dxc = size(Xc,2);           

% Set of Instruments
if d_iv == 1
    D = D(:,1); % only dist_port as IV
end

% Instrumental Variables
Z  = [D X];
dz = size(Z,2);    % Dimension of Z
    
% All Instruments -- Including constant term
Zc  = [Z ones(T,1)];
dzc = size(Zc,2);   % Dimension of Zc

% Vectors That Are Not Needed Anymore (Exogenous Regressors)
clear Alt Temp Rain Slope Soil  Mining Powerplant Pop Title
clear G_IR G50 G75 G100 G125 G150    


%% (B) ESTIMATE LAND-USE MODELS (LOGIT)


% I. Estimate Land-Use Regression Models: OLS and 2SLS

% (a) Logit OLS

% Coefficient
beta_ols = (Xc'*Xc)\(Xc'*lnY);

% Residuals
Uols = lnY - Xc*beta_ols;
  
% Robust Variance Matrix
Var_ols = inv(Xc'*Xc)*(Xc'*diag(Uols.^2)*Xc)*inv(Xc'*Xc);

% Standard Error
se_ols = sqrt(diag(Var_ols));

% t-statistic
t_ols = beta_ols./sqrt(diag(Var_ols));

% 95%-Confidence Interval
LB_ols = beta_ols - critical * sqrt(diag(Var_ols));
UB_ols = beta_ols + critical * sqrt(diag(Var_ols));


% (b) Logit 2SLS

% Projection Matrix
Pz = Zc*((Zc'*Zc)\Zc');

% Projection on X
Xc_hat = Pz*Xc;

% Coefficients
beta_iv = (Xc_hat'*Xc_hat)\(Xc_hat'*lnY);

% Residuals
Uiv = lnY - Xc*beta_iv;

% Robust Variance Matrix
Var_iv = inv(Xc_hat'*Xc_hat)*(Xc_hat'*diag(Uiv.^2)*Xc_hat)*inv(Xc_hat'*Xc_hat);

% Standard Error
se_iv = sqrt(diag(Var_iv));

% t-statistic
t_iv = beta_iv./sqrt(diag(Var_iv));

% 95%-Confidence Interval
LB_iv = beta_iv - critical * sqrt(diag(Var_iv));
UB_iv = beta_iv + critical * sqrt(diag(Var_iv));

% Overidentification Test
Jstat   = (Zc'*Uiv)' * inv(Zc'*diag(Uiv.^2)*Zc) * (Zc'*Uiv);
p_value = 1 - chi2cdf(Jstat,dzc-dxc); 


% III. Estimate Land-Use Quantile Regression Models: QR and IVQR

% (a) Logit QR
for i=1:t
    
    % Betas
    beta_qr_koencker = rq(Xc,lnY,tau(i));   
    
    % Store Logit QR betas
    beta_qr(:,i) = beta_qr_koencker; 
    
    % Standard Errors
    b_qr = vcqr(beta_qr(:,i),lnY,[TC X],tau(i));        
        
    % Store Logit QR standard errors
    se_qr(:,i)   = b_qr(:,2);
    
    % t-statistic
    t_qr(:,i) = beta_qr(:,i)./se_qr(:,i);
    
    % Confidence Interval
    LB_qr(:,i) = beta_qr(:,i) - critical * se_qr(:,i);    % LB of CI
    UB_qr(:,i) = beta_qr(:,i) + critical * se_qr(:,i);    % UB of CI    
    
end

clear beta_qr_koencker b_qr


% (b) Logit IVQR (Chernozhokov and Hansen's IVQR) 
for i=1:t
    
    %disp('Quantile')
    %disp(i)
    
    % Beta
    beta_ch = inv_qr(lnY,TC,X,D,tau(i));   % Chernozhukov and Hansen's IVQR
    
    % Store Logit IVQR betas
    beta_ivqr(:,i) = beta_ch;
    
    % Standard Error
    [b_ivqr,vc]  = vciqr_oid(beta_ivqr(:,i),lnY,TC,X,D,tau(i));
    se_ivqr(:,i) = b_ivqr(:,2);
    
    % t-statistic
    t_ivqr(:,i) = beta_ivqr(:,i)./se_ivqr(:,i);
    
    % Confidence Interval
    LB_ivqr(:,i) = beta_ivqr(:,i) - critical * se_ivqr(:,i);    % LB of CI
    UB_ivqr(:,i) = beta_ivqr(:,i) + critical * se_ivqr(:,i);    % UB of CI
    
    % IF estimate SPQIV Model, adjust t-statistic for ratio [-beta./beta(1)]
    % using the delta method for an appropriate comparison.        
    if est_spqiv == 1
        var_r = zeros(dxc-1,1);
        t_r   = zeros(dxc-1,1);
        
        for j=2:dxc-1
            % Variance of ratio from Delta Method
            var_r(j,1) = vc(1,1) * (beta_ivqr(j,i)^2/beta_ivqr(1,i)^4)...
                - 2 * vc(1,j) * (beta_ivqr(j,i)/beta_ivqr(1,i)^3)...
                + vc(j,j) * (1/beta_ivqr(1,i)^2);
            % Standard Error
            se_r(j,1) = sqrt(var_r(j,1));
            % t-statistic for ratio
            t_r(j,1) = (-beta_ivqr(j,i)/beta_ivqr(1,i))/sqrt(var_r(j,1));
        end
        % Store Standard Error for ratios
        se_ivqr_ratio(:,i) = se_r;
        % Store t-statistic for ratios
        t_ivqr_ratio(:,i) = t_r;        
        clear vc var_r t_r se_r
    end        
       
end

clear beta_ch b_ivqr


% IV. Marginal Effects at Sample Average (for 2SLS and IVQR)

% 2SLS
% a. Marginal Effects for TC
% Share of Deforestation at Covariates Sample Mean and Quantile tau
y_iv = exp(mean(Xc) * beta_iv) ./ (1 + exp(mean(Xc) * beta_iv));

% Derivative
dy_iv_tc = beta_iv(1) * y_iv * (1 - y_iv);

% b. Marginal Effects of Covariates at the Median 
dy_iv_cov = beta_iv * y_iv * (1 - y_iv);
    
% IVQR    
for i=1:ts_length
    
    % a. Marginal Effects for TC
    % Share of Deforestation at Covariates Sample Mean and Quantile tau
    y_ivqr(:,i) = exp(mean(Xc) * beta_ivqr(:,ts(i))) ./ (1 + exp(mean(Xc) * beta_ivqr(:,ts(i))));
    
    % Derivative
    dy_ivqr_tc(:,i) = beta_ivqr(1,ts(i)) * y_ivqr(:,i) * (1 - y_ivqr(:,i));
    
    % b. Marginal Effects of Covariates at the Median 
    dy_ivqr_cov(:,i) = beta_ivqr(:,ts(i)) * y_ivqr(:,i) * (1 - y_ivqr(:,i));
       
end


% V. Rearrange Quantiles to Avoid Quantile Crossing

% Vector with the municialities' unobservables
Uqr   = zeros(T,1);
Uivqr = zeros(T,1);

for j=1:T
    
    % (a) QR
    % Compute the quantile treatment response function for observation j
    q = Xc(j,:) * beta_qr;
    
    % Rearrange quantiles for obs j (to avoid quantile crossing)
    [q, ~] = sort(q);
    m = 1;
    while (lnY(j) > q(m) && m < size(q,2))    
        m = m + 1;
    end
    
    % Unobservable (quantile level) for observation j
    Uqr(j,1) = m;
    
    % OBS:  If lnY is between two estimated quantiles, choose the
    %       larger quantile.
    
    % (b) IVQR
    % Compute the quantile treatment response function for observation j
    q = Xc(j,:) * beta_ivqr;
    
    % Rearrange quantiles (to avoid quantile crossing)
    [q , ~] = sort(q);
    m = 1;
    while (lnY(j) > q(m) && m < size(q,2))    
        m = m + 1;
    end
    
    % Unobservable (quantile level) for observation j
    Uivqr(j,1) = m;
    
    % OBS:  If lnYs is between two estimated quantiles, choose the 
    %       larger quantile.
            
end



%% (C) SAVE ESTIMATION RESULTS


% I. Save the Data and the Estimation Results

   
if farm == 1
    
    % (a) Estimated Coefficients: Land-Use Regression
    beta_ols_farm_1      = [beta_ols se_ols t_ols LB_ols UB_ols];
    beta_iv_farm_1       = [beta_iv  se_iv  t_iv  LB_iv  UB_iv];
    beta_qr_farm_1       = [beta_qr(:,ts)   se_qr(:,ts)   t_qr(:,ts)   LB_qr(:,ts)   UB_qr(:,ts)];
    beta_ivqr_farm_1     = [beta_ivqr(:,ts) se_ivqr(:,ts) t_ivqr(:,ts) LB_ivqr(:,ts) UB_ivqr(:,ts)];
    pvalue_overid_farm_1 = p_value;
    obs_farm_1           = T;
    margeff_farm_1       = [100*[y_iv; dy_iv_tc * 10],...
                            100*[y_ivqr(1); dy_ivqr_tc(1) * 10],...
                            100*[y_ivqr(2); dy_ivqr_tc(2) * 10],...
                            100*[y_ivqr(3); dy_ivqr_tc(3) * 10],...
                            100*[y_ivqr(4); dy_ivqr_tc(4) * 10],...
                            100*[y_ivqr(5); dy_ivqr_tc(5) * 10]];

    % Save All
    save beta_ols_farm_1      beta_ols_farm_1      -ascii
    save beta_iv_farm_1       beta_iv_farm_1       -ascii
    save beta_qr_farm_1       beta_qr_farm_1       -ascii
    save beta_ivqr_farm_1     beta_ivqr_farm_1     -ascii
    save pvalue_overid_farm_1 pvalue_overid_farm_1 -ascii
    save obs_farm_1           obs_farm_1           -ascii
    save margeff_farm_1       margeff_farm_1       -ascii
    
    % (b) If run SPQIV, save t_ratios
    if est_spqiv == 1
        t_ratio_farm1 = t_ivqr_ratio(:,50);
        save t_ratio_farm1  t_ratio_farm1   -ascii
    end
    
    % (c) Unobservables and Coefficients  
    Data_farm_1 = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)'];
    save Data_farm_1   Data_farm_1   -ascii
        
    
elseif farm == 2
    
    % (a) Estimated Coefficients: Land-Use Regression
    beta_ols_farm_2      = [beta_ols se_ols t_ols LB_ols UB_ols];
    beta_iv_farm_2       = [beta_iv  se_iv  t_iv  LB_iv  UB_iv];
    beta_qr_farm_2       = [beta_qr(:,ts)   se_qr(:,ts)   t_qr(:,ts)   LB_qr(:,ts)   UB_qr(:,ts)];
    beta_ivqr_farm_2     = [beta_ivqr(:,ts) se_ivqr(:,ts) t_ivqr(:,ts) LB_ivqr(:,ts) UB_ivqr(:,ts)];
    pvalue_overid_farm_2 = p_value;
    obs_farm_2           = T;
    margeff_farm_2       = [100*[y_iv; dy_iv_tc * 10],...
                            100*[y_ivqr(1); dy_ivqr_tc(1) * 10],...
                            100*[y_ivqr(2); dy_ivqr_tc(2) * 10],...
                            100*[y_ivqr(3); dy_ivqr_tc(3) * 10],...
                            100*[y_ivqr(4); dy_ivqr_tc(4) * 10],...
                            100*[y_ivqr(5); dy_ivqr_tc(5) * 10]];
    
    % Save All
    save beta_ols_farm_2      beta_ols_farm_2      -ascii
    save beta_iv_farm_2       beta_iv_farm_2       -ascii
    save beta_qr_farm_2       beta_qr_farm_2       -ascii
    save beta_ivqr_farm_2     beta_ivqr_farm_2     -ascii
    save pvalue_overid_farm_2 pvalue_overid_farm_2 -ascii
    save obs_farm_2           obs_farm_2           -ascii
    save margeff_farm_2       margeff_farm_2       -ascii
    
    % (b) If run SPQIV, save t_ratios
    if est_spqiv == 1
        t_ratio_farm2 = t_ivqr_ratio(:,50);
        save t_ratio_farm2  t_ratio_farm2   -ascii
    end
    
    % (c) Unobservables and Coefficients
    Data_farm_2 = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)'];
    save Data_farm_2  Data_farm_2   -ascii
    
elseif farm == 3
    
    % (a) Estimated Coefficients: Land-Use Regression
    beta_ols_farm_3      = [beta_ols se_ols t_ols LB_ols UB_ols];
    beta_iv_farm_3       = [beta_iv  se_iv  t_iv  LB_iv  UB_iv];
    beta_qr_farm_3       = [beta_qr(:,ts)   se_qr(:,ts)   t_qr(:,ts)   LB_qr(:,ts)   UB_qr(:,ts)];
    beta_ivqr_farm_3     = [beta_ivqr(:,ts) se_ivqr(:,ts) t_ivqr(:,ts) LB_ivqr(:,ts) UB_ivqr(:,ts)];
    pvalue_overid_farm_3 = p_value;
    obs_farm_3           = T;
    margeff_farm_3       = [100*[y_iv; dy_iv_tc * 10],...
                            100*[y_ivqr(1); dy_ivqr_tc(1) * 10],...
                            100*[y_ivqr(2); dy_ivqr_tc(2) * 10],...
                            100*[y_ivqr(3); dy_ivqr_tc(3) * 10],...
                            100*[y_ivqr(4); dy_ivqr_tc(4) * 10],...
                            100*[y_ivqr(5); dy_ivqr_tc(5) * 10]];
    
    % Save All
    save beta_ols_farm_3      beta_ols_farm_3      -ascii
    save beta_iv_farm_3       beta_iv_farm_3       -ascii
    save beta_qr_farm_3       beta_qr_farm_3       -ascii
    save beta_ivqr_farm_3     beta_ivqr_farm_3     -ascii
    save pvalue_overid_farm_3 pvalue_overid_farm_3 -ascii
    save obs_farm_3           obs_farm_3           -ascii
    save margeff_farm_3       margeff_farm_3       -ascii
    
    % (b) If run SPQIV, save t_ratios
    if est_spqiv == 1
        t_ratio_farm3 = t_ivqr_ratio(:,50);
        save t_ratio_farm3  t_ratio_farm3   -ascii
    end
    
    % (c) Unobservables and Coefficients
    Data_farm_3 = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)'];
    save Data_farm_3   Data_farm_3   -ascii
    
elseif farm == 4
    
    % (a) Estimated Coefficients: Land-Use Regression
    beta_ols_farm_4      = [beta_ols se_ols t_ols LB_ols UB_ols];
    beta_iv_farm_4       = [beta_iv  se_iv  t_iv  LB_iv  UB_iv];
    beta_qr_farm_4       = [beta_qr(:,ts)   se_qr(:,ts)   t_qr(:,ts)   LB_qr(:,ts)   UB_qr(:,ts)];
    beta_ivqr_farm_4     = [beta_ivqr(:,ts) se_ivqr(:,ts) t_ivqr(:,ts) LB_ivqr(:,ts) UB_ivqr(:,ts)];
    pvalue_overid_farm_4 = p_value;
    obs_farm_4           = T;
    margeff_farm_4       = [100*[y_iv; dy_iv_tc * 10],...
                            100*[y_ivqr(1); dy_ivqr_tc(1) * 10],...
                            100*[y_ivqr(2); dy_ivqr_tc(2) * 10],...
                            100*[y_ivqr(3); dy_ivqr_tc(3) * 10],...
                            100*[y_ivqr(4); dy_ivqr_tc(4) * 10],...
                            100*[y_ivqr(5); dy_ivqr_tc(5) * 10]];
    
    % Save All
    save beta_ols_farm_4      beta_ols_farm_4      -ascii
    save beta_iv_farm_4       beta_iv_farm_4       -ascii
    save beta_qr_farm_4       beta_qr_farm_4       -ascii
    save beta_ivqr_farm_4     beta_ivqr_farm_4     -ascii
    save pvalue_overid_farm_4 pvalue_overid_farm_4 -ascii
    save obs_farm_4           obs_farm_4           -ascii
    save margeff_farm_4       margeff_farm_4       -ascii
    
    % (b) If run SPQIV, save t_ratios
    if est_spqiv == 1
        t_ratio_farm4 = t_ivqr_ratio(:,50);
        save t_ratio_farm4  t_ratio_farm4   -ascii
    end
    
    % (c) Unobservables and Coefficients
    Data_farm_4 = [Code Y Area Qa1 Qa2 Uivqr beta_ivqr(1,Uivqr)'];
    save Data_farm_4    Data_farm_4   -ascii
    
end



%% (D) DELETE IF NOT NEEDED

% Estimates
clear Var_ols Var_iv se_ols se_iv se_qr se_ivqr  t_ols t_iv t_qr t_ivqr
clear LB_ols LB_iv LB_qr LB_ivqr  UB_ols UB_iv UB_qr UB_ivqr
clear Uiv_hat R2 J critical_J p_value

% Extra Stuff
clear i index j m p q 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
